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Abstract 

We present a compact discontinuous Galerkin (CDG) method for an 
elliptic model problem. The problem is first cast as a system of first order 

equations by introducing the gradient of the primal unknown, or flux, as 
an additional variable. A standard discontinuous Galerkin (DG) method is 
then applied to the resulting system of equations. The numerical interele- 
ment fluxes are such that the equations for the additional variable can be 
eliminated at the element level, thus resulting in a global system that in- 
volves only the original unknown variable. The proposed method is closely 
related to the local discontinuous Galerkin (LDG) method [B. Cockburn 
and C.-W. Shu, SIAM J. Numer. Anal, 35 (1998), pp. 2440-2463], but, 
unlike the LDG method, the sparsity pattern of the CDG method in- 
volves only nearest neighbors. Also, unlike the LDG method, the CDG 
method works without stabilization for an arbitrary orientation of the el- 
ement interfaces. The computation of the numerical interface fluxes for 
the CDG method is slightly more involved than for the LDG method, 
but this additional complication is clearly offset by increased compact- 
ness and flexibility. Compared to the BR2 [F. Bassi and S. Rebay, J. 
Comput. Phys., 131 (1997), pp. 267-279] and IP [J. Douglas, Jr., and 
T. Dupont, in Computing Methods in Applied Sciences (Second Internal. 
Sympos., Versailles, 1975), Lecture Notes in Phys. 58, Springer, Berlin, 
1976, pp. 207-216] methods, which are known to be compact, the present 
method produces fewer nonzero elements in the matrix and is computa- 
tionally more efficient. 

1 Introduction 

Discontinuous Galerkin (DG) methods [11] have become the subject of con- 
siderable research over recent years due to their potential to overcome some 
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of the perceived shortcomings of the more established discretization methods. 
For convection problems, DG methods produce stable discretizations without 
the need for cumbersome stabilization strategies. They work well on arbitrary 
meshes and allow for different orders of approximation to be used on different 
elements in a very straightforward manner. Clearly, this flexibility comes at the 
expense of duplicating the degrees of freedom at the element boundary inter- 
faces. This is a serious drawback when low order polynomial approximations 
are used, but it is less important for high order interpolations. DG methods 
appear to be ideally suited for applications involving wave propagation phenom- 
ena, where low dispersion and high accuracy are required, such as aeroacoustics 
or electromagnetics. 

While DG methods seem to be well suited for the discretization of first order 
hyperbolic problems, their extension to elliptic problems is far less obvious. A 
number of extensions to deal with the elliptic problem have been proposed and 
analyzed under a unified framework in [1]. Also, a comparison of the perfor- 
mance of various schemes from a practical perspective is presented in [6] . Among 
the various alternatives, the local discontinuous Galerkin (LDG) method [10] has 
emerged as one of the most popular choices. The LDG method appears to be one 
of the most accurate and stable schemes among those tested. In addition, the 
LDG method is easy to implement for complex convective-diffusive systems and 
can be generalized to handle equations involving higher order derivatives [16]. 
In the LDG method, the original equation involving second order derivatives is 
cast as a system of first order equations by introducing additional variables for 
the solution gradient, or fiux. The resulting system is then discretized using 
a standard DG approach. By appropriately choosing the interelement fluxes, 
the additional variable can be eliminated locally. Thus, a stable discretization 
that involves only the original unknown variable is obtained. Unfortunately, 
when the LDG method is used in miiltiple dimensions, the discretization gener- 
ated has the undesirable feature that the degrees of freedom in one element are 
connected, not only to those in the neighboring elements, but also to those in 
some elements neighboring the immediate neighbors. For applications employ- 
ing explicit or iterative solution techniques, this is usually not a problem, but 
for applications where the matrix needs to be formed, this represents a severe 
disadvantage. 

Two alternative formulations for the treatment of the second order deriva- 
tives are the symmetric interior penalty (IP) method [12] and the BR2 method 
proposed in [3]. In these methods, the original form of the equation involving 
second derivatives is discretized directly, and stabilization is added explicitly in 
a suHicient amount to render the method stable. Although somewhat simpler, 
the IP method appears to be less popular than the BR2 method. This is prob- 
ably because of the requirement of a penalty parameter that depends on both 
the mesh and the approximation order. Both these methods have the advantage 
that they are compact in the sense that only the degrees of freedom belonging 
to neighboring elements are connected in the discretization. When suitable pe- 
nalization is employed these approaches are competitive with the LDG scheme 
in terms of accuracy. Thus, these schemes are an attractive alternative to the 
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LDG scheme when an impHcit solution of the discretized system is required. 

For many applications of interest involving convective-diffusive systems, such 
as the Navier-Stokes equations at high Reynolds numbers, the time and length 
scales are such that implicit discretization turns out to be a requirement. In this 
paper, we develop a variation of the LDG method, the compact discontinuous 
Galerkin (CDG) method. The main motivation for developing this new scheme 
is to eliminate the distant connections between nonneighboring elements which 
arise when the LDG scheme is used in multiple dimensions. We note that in 
the one-dimensional case the CDG and LDG schemes are identical, but in the 
multidimensional case they differ in the approximation to the solution gradient 
at the interface between neighboring elements. This seemingly minor difference 
results in a scheme that appears to inherit all the attractive features of the LDG 
method and is compact. In addition, numerical experiments indicate that the 
CDG scheme is slightly more stable than the LDG method and is less sensitive 
to the element and/or interface orientation. In particular, when the stabilization 
constant is set to zero, the CDG scheme is stable in situations where the LDG 
method is unstable. It is well known that, without explicit stabilization, the 
LDG scheme is stable only when the orientation of element interfaces satisfies 
a certain condition [15]. 

Since the CDG scheme is compact, it produces a sparser connectivity matrix 
than the LDG scheme, meaning lower storage requirements and higher computa- 
tional performance. Thus, the slight additional increase in complexity involved 
in the numerical flux evaluation is more than offset by the increased efficiency 
benefits. Compared to the IP and BR2 methods, the CDG scheme is computa- 
tionally simpler, generates a sparser matrix with a smaller number of nonzero 
elements when using a nodal basis, and appears to produce slightly more accu- 
rate results than the BR2 method in the numerical tests performed. Given the 
similarities between the BR2 and IP methods, we have considered only the BR2 
method in our numerical comparisons. 

The remainder of the paper is organized as follows. In section 2, we intro- 
duce our model second order elliptic problem. Next, we describe the LDG dis- 
cretization method and adopt the framework introduced in [1] to write the LDG 
algorithm in the so-called primal form. This form, involving only the original 
problem variable, highlights the symmetry of the scheme as well as the sparsity 
pattern. In section 3, we present the CDG method. The CDG method is then 
written in primal form so that it can be easily compared with the LDG method. 
Like the LDG method, the CDG method is shown to be symmetric, conserva- 
tive, and adjoint consistent. It turns out that the CDG and LDG schemes are 
so closely related that the error estimate presented in [1] for the LDG method 
is essentially applicable to the CDG method without changes. In section 4, we 
compare the LDG and CDG schemes using the test problem presented in [15]. 
The increased stability of the CDG scheme, for arbitrary interface ordering, is 
shown numerically by calculating the size of the null-space for the model test 
problem. Practical implementation and efficiency issues such as sparsity pat- 
terns and storage requirements for the LDG, BR2, and CDG schemes, in the 
more general d-dimensional setting, are addressed in section 5. Finally, we con- 
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elude in seetion 6 with some numerical results aimed at comparing the accuracy 
and conditioning of the LDG, BR2, and CDG schemes. 



2 Discontinuous Galerkin formulation 

2.1 Problem definition 

The proposed method will be described for the model Poisson problem 
-V • (kVu) = / in n, 

u = qd on dQ.D, (^i-^ 

where Q is a bounded domain in with boundary 90 = 50 £> U 50 jv and 
d = 1,2, or 3 is the dimension. Here, f{x) is a given function in L?{Q), and 
n{x) e L°°{Q) is positive. Further, we assume that the length of dQ.D is not 
zero. 



2.2 DG Formulation for elliptic problems 

In order to develop a DG method, we rewrite the above problem (1) as a first 
order system of equations 
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where n is the outward unit normal to the boundary of O. 

Next, we introduce the broken spaces V{Tfi) and '^{Th) associated with the 
triangulation 7^ = {K} of O. In particular, V{Tf^) and ^{Th) denote the spaces 
of functions whose restriction to each element K belongs to the Sobolev spaces 
H^iK) and [H\K)]'^. That is, 

V = {v G L'^ifl) \ v\k & H\K) yX GTh}, (3) 

s = {Te[L\Q)f\T\Ke[H\K)fyKeTh}. (4) 

In addition, we introduce the finite element subspaces Vh CV and S/j c S as 

Vh = {ve L\n) I v\k e Vp{K) yK e (5) 
Eh = {re [L\n)]'' I r\K e [Vp{K)f e %}, (6) 

where 'Pp{K) is the space of polynomial functions of degree at most p>lonK. 
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Following [10], wc consider DG formulations of the form: find Uh S Vh and 
(Th G S?t such that for all G 7^ we have 



/ (Th ■ T dx = ~ UhV ■ (kt) dx+ UKT -nds Vr e ['Pp(iV')]'^, (7) 
Jk Jk JdK 

/ (Jh-Vvdx = - fvdx+ a-nvds \/veVp{K). (8) 
Jk Jk JdK 

Here, the numerical fluxes & and u are approximations to cr = kVu and to 
u, respectively, on the boundary of the element K. The DG formulation is 
complete once we specify the numerical fluxes & and u in terms of (Th and Uh 
and the boundary conditions. 

Expressions (7) and (8) apply to each clement separately. In order to write 
expressions which are applicable over the whole domain, we require some addi- 
tional notation. Here, we closely follow the notation used in [1]. 

Gonsider two adjacent elements and of the triangulation 7;,, and 
denote by e = dK^ r\dK~ their common face. Further, assume that denote 
the unit normals to dK^, respectively, at any point on the face e. Similarly, let 
(t*, denote the traces on e of functions (r, v) E E;i x Vh which are smooth 
in the interior of elements . The average and jump operators are given as 

{t} = (t+ + t-)/2, {v} = {v+ + v-)/2, 

[t] = r+ • n+ + T~ ■ n~, [v] = + v~n~ . 

Note that, according to this definition, the jump of a scalar quantity is a vector, 
but the jump of a vector quantity becomes a scalar. 

Now, by summing (7) and (8) over all elements and considering only con- 
servative schemes for which the numerical fiuxes u and d- on a given face are 
unique, we obtain the following global expressions: find Uh & Vh and (Th € ^h 
such that 

/ (Th ■ T dx = — UhVh ■ (kt) dx + u[kt] ds + ukt ■ nds Vr e E/j, 
Jn Jn Jsi Jan 

(9) 



(10) 



/ (Th-^hvdx = / fvdx+ / &-[v]ds+ / va-nds e V/j, 
Jn Jn Jsi Jan 

where £i denotes the union of all the interior faces in the triangulation 7^ . Also, 
Vft, denotes the broken gradient operator. That is, V/jU and V/j, • r are functions 
whose restriction to K is equal to Vu and V • r, respectively. 

For later use, we note that, if we use the integration by parts formula, 

— / vVh-'i'dx = I T-S/hvdx— / ([«] • {t} + {w}[r]) ds— / VT-nds, (11) 
Jn Jn Jsi Jan 
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which is vaUd for all r e [H^{Th)]'^ and v G H^{Th), we can write (9) as 

/ CTh-Tdx= / T ■ {kV hUh) dx - / {[uh] ■ {kt} - {u - Uh}[KT]) ds 
Jn Jn JSi 

(12) 

+ / {u — UhjKT ■ nds VrGTih- 

Jon 

2.3 The LDG method 

Since our method is c;losely related to the LDG method presented in [10], we 
start with a description of the LDG algorithm. For the LDG method, the 
numerical interelement fluxes {a, u) are given by 

a = W}-CiiK] + Ci2K], (13) 
U = {Uh} - Ci2 ■ [uh] (14) 
for the interior faces, and 

a-^(Th- Cii{uh- gD)n, u = 9d on dVlo, 

& = gNTi, u = Uh on dQ,N, 



(15) 



for the boundary faces. Here, Cn is a positive constant and C12 is a vector 
which is determined for each interior face according to 

C^2=\{S^ln+ + S^ln-), (16) 



where e {0, 1} is a switch which is defined for each element face. That is, 
denotes the switch associated with element on the face that element 
shares with element K~ . The switches always satisfy that 



S^l+S^l=l (17) 

but are otherwise arbitrary. We note that, although the form (16) is not the 
most general form for C12 presented in [10], other choices lead to wider stencils 
in the final discrete equations. We also point out that the choice of element face 
switches has an effect on the final form of the discrete equations. 

2.3.1 Primal form of the LDG algorithm 

In order to derive the primal form of the LDG algorithm, we first particularize 
(12) for the fluxes given by (14), 



/ ah-Tdx= / T • {nVhUh) dx - {[uh] ■ {kt} + C12 ■ [wh][«T]) ds 
Jn Jn JSi 

+ {go- Uh)KT ■ nds Vr e E^. 

Jdnn 



(18) 
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To obtain an expression for cr^ as a function Uh, we follow [1] and introduce the 
lifting operators r : [L^(fi)]'' ^ ^h, I ■ L^{£i) ^ ^h, and td : L^idno) ^h- 

/ r{(i))-Tdx = - (j)-{T}ds VreEh, 
Jn J Si 

I l{q)-Tdx = - q[T]ds Vr e S,,, (19) 

Jn JEi 

/ rD{q)-Tdx = — QT-nds Mt &T,h- 
Jn Jann 

Thus, we can write (18) as 

/ {a-h - nVhUh - K,r{[uh]) - kI{Ci2 ■ [uh]) + uroigo - uu)) • rdx = 

(20) 

Vr e S^. 

Therefore, we have 

CTh = KVhUh + ^h, (21) 

where e S/, is 

d-h = Kr{[uh]) + kI{Ci2 ■ [uh]) - Kroigo - Uh). (22) 

Thus, we see that that cr^ is equal to KVhUh plus an additional perturbation 
term which is forced by [uh], Ci2-[uh], and go — Uh- Also, note that roigD — Uh) 
is nonzero only on the elements that have a face on the Dirichlct boundary. In 
writing expressions (21) and (22), we have assumed that V/jV?i C S/^, which is 
certainly the case if equal order polynomial interpolants are used for 14 and S/j. 
Setting T = VhV in (18), we can rewrite (10) as 

/ V hV ■ {k\/ hUh) dx - / {[uh] ■ {k"^ hv} + Ci2 ■ [uh][K\/ hv]) ds 
Jn JSi 

+ / {go - Uh)nVhV ■ nds 
Jdno 

= fvdx + a--[v]ds+ va-nds VveVh- (23) 
Jn Jsi Jon 

Making use of (13), (15), (21), and (22), the terms involving <t in the above 
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equation can be written as 

/ &-[v]ds^ i{KVhUh} + Ci2[KyhUh]) -^(13+ / i{a} + Ci2[^])-[v]ds 

= j {{nWhUh) + Ci2[KWhUh]) ■ [v] ds 

- [ + KC12 ■ [v])) ■ {r{[uh]) + l{Ci2 ■ [uh]) + TDiuh)) dx 

Jn 

[ '^{r{[v]) + l{Ci2 ■ [v])) ■ rn{9D) dx - [ Cn[u]-[v]ds 
Jn JSi 



+ 
and 



/ va nds= vcrh-nds— / Ciivuhds+ / CiivgDds+ / vqn ds 

Jdn JdQn JdQD Jdno J dn^ 

= I VKVhUh ■ nds + / Kv{r{[uh]) 

JdQn JdQn 



+ /(C12 • [uh]) +rD{uh)) ■ nds 



- / KvroigD) ■ nds + / C-iiv{gD - Uh) ds + / vgnds 

J dQo J dflo J OQm 

/ VKVhUh ■ nds — / uroiv) 
Jdno Jn 



ldQ.D 

{^{[uh]) + l{Ci2 ■ [uh]) + roiuh)) dx 



- / nvroigo) ■ nds + / Ciiv{gD - Uh) ds + / vg^ds. 

JdQn J dnn JoQm 

Therefore, we can rewrite (23) as 

B^^«K,z;) = L^^«(^;) yveV^, (24) 
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where the bilinear form Bf^^*^ : 14 x V/j — > M is given by 

B^^^'{u,v)= / y hv ■ {kS/ hu) dx - / {[u] ■ {kS/ hv} + {nV hu} ■ [v]) ds 



- j (Ci2 • [u] [nVhv] + [kVam]Ci2 ■[v])ds + j Cn [u] ■ [v] ds 

n{r{[u]) + /(Ci2 • [u]) + rD{u)) ■ {r{[v]) + l{Ci2 ■ [v]) + W) dx 

^ / {k^ hU ■ nv + ukS/ hV ■ n)ds + / Cnuvds (25) 
and the Unear form Lj^^'^ : M is given by 

Lh^^iv) = [ fvdx- f gninVhV + r{[v]) + Z(Ci2 • [v])) ■ nds 

- / KvrD{gD)-nds+ j Ci\gDvds+ j vgnds Vt; e V/,. 

(26) 

It is straightforward to verify that the biUnear form (25) is symmetric, i.e., 

Bh{u,v) ~ Bh{v,u). Also, the conservative form of the numerical fluxes, (13) 
and (14), guarantees that the LDG scheme is conservative and adjoint consistent 
[!]• 

Unfortunately, when the scheme is implemented in multidimcnsions on gen- 
eral triangular/tetrahedral meshes, the resulting discretization is not compact 
in the sense that the equation corresponding to a given degree of freedom may 
involve degrees of freedom that belong to elements which arc not immediate 
neighbors. It turns out that these additional connections are due to the volume 
term in (25) which involves products of the lifting functions. Although the con- 
nectivity pattern between elements depends on the choice of face switches in 
(16), it is well known [15] that in multidimcnsions this problem cannot be reme- 
died by a more careful choice of the face switches (16). This noncompactness of 
the LDG scheme occurs also for quadrilateral/hexahedral discretizations. 



3 The CDG algorithm 

The CDG algorithm is designed to be compact and, at the same time, inherit all 
the attractive properties of the LDG algorithm. To start with, we decompose 
the lifting operators introduced in (19) into facewise contributions. Thus, we 
consider for all e G £i, r"^ : [L^{e)]'^ — > S^, : i^(e) — > T,h and for each 
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e e dflo, I'D '■ -t'^(e) — > S/i, defined as 

[ r''{(f))-Tdx = -[^-{rjds Vr e S^, 
in Je 

[ r{q)-Tdx = -[q[T]ds VreSfc, (27) 

in Je 

/ fDil)'''''^^ = — qr-nds Vr G S/j. 

Jn 

Clearly, we will have, for all (p e [//^(f,)]'' and all q € L^{£i), 

r{ct>) = J2r%4>), l{q) = J2nQ), Mq) = J2 "-Dil)- (28) 

Now, we can define the CDG method. The numerical interelement fltixes 
((T, u) for the CDG method are given by 

&^{cTl}-Cn[uH] + Ci2[cTl], (29) 

u = {uh} - Ci2 ■ [uh] (30) 



(31) 



for the interior faces, and 

a = al - Cn{uh - gD)n, u = gD on dfln, 

& = gpfn, u = Uh on dQ,N, 

for the boundary faces. Here, cr^ is given as 

al = KVhUh + ^l, (32) 

where 

= Kr''{[uh]) + Kr{Ci2 ■ [uh]) - nrf^igo - Uh) ■ (33) 

We note that the numerical flux, u, is chosen as in the LDG method. Therefore, 
(18) and (21)-(23) still apply for the CDG method, and the only difference 
between the LDG and CDG methods is in the evaluation of the terms involving 
a in (23), which in the CDG case is done according to (29) and (31). Also, the 
coefficients C12 are given by expressions (16) and (17). 

In order to compute the CDG numerical flux it on a given face e, we need 
to evaluate first a stress field af^ associated with this face. This evaluation, 
however, can be carried out efficiently due to the localized support of In 
particular, we note that when e € OUn, then o-£ = 0. When e € dil,]j, we 
have — Kr'jj^gu — Uh), which has only a nonzero support on the clement 
neighboring face e. Finally, when e € Si, then d-^ = Kr^{[uh]) + nl'^iCri ■ [uh])- 
In this case, a'f^ is nonzero only in one of the elements neighboring face e. The 
element in which ct£ is nonzero is determined by the choice of switches for that 
face. In particular, using (16) and (27), it can be easily shown that if 5^+ = 1 
and S^t. = 0, then = on K~ . Similarly, we will have ct^ = on when 
= and S§1 = 1. 
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3.1 Primal form of the CDG algorithm 

In order to obtain the primal form of the CDG method, we proceed as before 
and start from (23). In this case, the terms involving a become 



/ & ■[v]ds= / a -[v] ds 

= [ {{KVhUh} + Ci2[KVhUh]) -Mds+Y, [ + Ci2[*1) • [v] ds 

- [ Cn[u]-[v]ds 
= / {{K^hUh} + Ci2[KVhUh]) ■ [v] ds 

-Y. I <r''{['^]) + nCi2^v]))-{r\[uh]) + l'{C^2-[uh])+rUuh))dx 

+ V [ n{r-{[v]) + r{C^2-[v\))-r'i,{gD)dx- [ Cn[u]-[v]ds 



and 



/ v & ■ nds = ^ / v & ■ nds 

= / vafi^-nds— I Cuvuhds+ I Cwygods + I vgNds 

= VKVhUh-nds+ ^ Kv{r''{[uh])+l''{Ci2-[uh]) + r'iy{uh))-nds 

- ^ Kvr'jjigD) ■ nds - / Ciivuhds+ / CiivgDds+ / vgN ds 

= / VK\/hUh-nds- ^ / Kr^j^iv)- {r%[uh])+l''{Ci2-[uh]) + r'}j{uh))dx 
Jano eeaQo ^ 

~ X/ / i^'"1'd{9d) ■ nds - / CiiVUhds+ / CiivgDds+ / vgN ds. 

Thus, for the CDG scheme, (23) can be written as 

B^^'^iuh, v) = L^^^(^;) V^; G V^, (34) 
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where the bilinear form B^^*^ : V/j x V/i — > K is given by 

B^^°{u, v)= Vhv ■ {K'Vhu) dx- {[u] ■ {kWhv} + {kWhu} ■ [v]) ds 

JQ. JSi 

- / (Ci2 • M[kV^w] + [KVfcu]Ci2 • [v]) ds 

JSi 

+ E / «(r^(M) + «'(Ci2-M) + r!,(«))-(r^(M)+r(Ci2-M) + r|,(t;))rfa; 

— / (kV/(U • nv + unVhV ■ n)ds + / Cii[u] • [v] ds + Cnuv ds 

(35) 

and the linear form i*^^*^ : 14 — > M is given by 
Lh^^{v)= / fvdx- / goK^hV-nds 

JQ JdQo 

- I nvroigo) ■ nds+ / CiigDvds+ / vgNds WveVh. 

J OQd ^ OQd ^ OQm 

(36) 

The CDG method is symmetric, i.e., B^^'^{u,v) = B^^'^{v,u), and re- 
tains all the attractive properties of the LDG algorithm such as consistency and 
adjoint consistency. 

3.2 Error estimates 

We observe that the only difference between the LDG and CDG schemes is the 
stabilizing term involving the products of the lifting functions. In the LDG 
scheme, we have 

/ K{r{[u]) + l{Cu ■ [u]) + roiu)) ■ {r{[v]) + /(C12 • [v]) + roiv)) dx, 

JQ 

EE/ '^(^'(M) + • M) + r|,(n)) • {rf{[v]) + lf{C^2 ■ [v]) + r{,{v)) dx, 

eeSi feSi ''^ 

(37) 

whereas in the CDG scheme, we have 

E / <r%[u])+l%C,2-[u])+rUu))-{r%[v]) + l%C,2-[v]) + rh{v))dx, 
JQ 



eeSi 



E E -^^Z / '^('^'(M) + «^(Ci2 • M) + rUu)) ■ {rf{[v]) + 1^ {C,, ■ [v]) + r{,{v)) dx, 

(38) 
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where 5ef is the Kroncckcr delta. Thus, we sec that the CDG scheme can be 
regarded as the LDG algorithm with some terms turned off. We also note that 
the turned-off terms in the LDG algorithm are indefinite and hence are not 
guaranteed to contribute to the method's stability. The effect of using lifting 
functions in the CDG method which are associated with individual faces is to 
eliminate connectivities between nonneighboring elements. We note that an 
analogous approach was adopted in [3, 5] to render the BR2 scheme compact. 

It turns out that the proofs of coercivity and boundedness for the LDG 
method presented in [1] can be used here without change. This leads to optimal 
a priori estimates for the CDG method, 

\\\u-Uh\\\<ChF\u\p+^,n (39) 

and 

\\u-Uh\\o,a<ChP+^\u\j,+i^n . (40) 
Here, the norm 1 1 1 • 1 1 1 is given by 

111^111'= E \<K+Y.\Mv\)\\ln+ E II^^MIIo,Q- (41) 

The above estimates require that the stabilization parameter Cn in (29) 
is taken to be of order 0(h~^), where h is the characteristic mesh size (see 
also [7]). We note that for Cn of order 0(1), only suboptimal convergence is 
demonstrated, but in practical computations, optimal results are also observed. 
We also point out that for general discretizations, the piecewise constant approx- 
imation p = does not lead to a consistent discretization. This is in common 
with other DC schemes such as the LDG or the BR2. 



4 Stabilization 

The above a priori error estimates are applicable to both the CDG and LDG 
algorithms. It turns out that, for the LDG algorithm, one can set Cn = for 

all the internal intcrfac;es, provided the switches in (16) are chosen following a 
simple rule. That is, if the switches for each simplex element K satisfy that 

Y^S^'<d+l, (42) 

eedK 

where d is the problem dimension, then the scheme shows no degradation in 
performance and becomes extremely simple. This result was proven in [8]. In 
this case, the numerical flux u on a given internal face is taken to be the value 
of Uh on one of the neighboring elements, while the numerical flux a is taken 
to be the value of cr^ on the other neighboring element. The element used to 
calculate either w or it is determined by the value of switches on that face. The 
rule (42) guarantees that, when calculating the numerical fluxes on each face, 
the value of the solution on each element will be used, at least once, to set u on 
the element boundary, and, at least once, to set a on the element boundary. 



13 



Clearly, there is plenty of flexibility in choosing appropriate values for switches 
which satisfy the rule (42); see [9], for instance. Thus, provided that the rule 
(42) is satisfied, the LDG scheme converges at the optimal rate without the 
need for explicit stabilization. 

4.1 Null-space dimension 

We have found that while the rule (42) is essential in ensuring that the solution 
is unique for the LDG method, this requirement is not necessary for the CDG 
method. That is, for the CDG method we are able to set Cn = for all the 
internal faces and use any combination of switches with the only constraint 
given by (17). 

In order to illustrate this point, we adopt the two-dimensional test problem 
presented in [15]. We consider a square domain with periodic boundary condi- 
tions imposed on all sides. We perform a regular subdivision into four squares 
and tlic^n subdivide each square into two triangles. We look at approximations 
ranging from p = 1 to p = 7 and nodal basis functions with equally spaced 
nodes. We discretize the Laplacian operator using the CDG and the LDG al- 
gorithms with the parameter Cn set to zero and calculate the dimension of the 
null-space of the resulting matrix. 

We consider two different switches for both the LDG and CDG algorithms. 
The so-called consistent switch satisfies (42) , and here it is chosen using a proce- 
dure analogous to that presented in [9, 15]. We also consider the natural switch, 
which is based on element numbering and sets 5^+ = if the element num- 
ber is less than the element number K~ , and to 1 otherwise. This switch 
was first introduced in [2] in the context of interior point methods for elliptic 
problems. 

Because of the periodic boundary conditions, any solution will be undeter- 
mined up to a constant, and as a consequence, we expect a singular matrix with 
a null-space of dimension one. The computed dimension of the null-space for 
the different schemes, polynomial order interpolations, and switches is presented 
in Table 1. We note that while the LDG scheme gives the desired null-space 
dimension of one when the consistent switch is employed, the null-space dimen- 
sion grows with increasing p, when the natural switch is employed. This same 
result was reported in [15]. On the other hand, the CDG scheme always gives 
the desired one-dimensional null-space for all p and for any switch choice. 

We note that the natural switch has some computational advantages when 
computing the ILU(O) factorization of the system matrix [14]. If = 

when < K~ , the lower triangular blocks in the matrix have only a few 
nonzero rows, and no additional fill-in is introduced during the factorization 
phase. On the other hand, for an arbitrary switch choice, some lower triangular 
blocks will have nonzero columns that will render the blocks completely full 
after factorization. This effect is described in more detail in [14], where the 
CDG method is used to discretize convective-diffusive systems which are solved 
using a preconditioned Krylov solver. 
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NuUspace dimension 



Polynomial order p 



1 2 3 4 5 6 7 



Consistent switch 



CDG 
LDG 
CDG 
LDG 



1111111 
1111111 
1111111 
3 4 5 6 7 8 9 



Natural switch 



Table 1: Nullspace dimensions for the CDG/LDG schemes using the two differ- 
ent switches. The problem is expected to have a one-dimensional null-space, but 
with the (inconsistent) natural switch the LDG scheme gives spurious modes 
and a null-space that grows with p. 



5 Implementation 

Since the main motivation for developing the CDG algorithm is to obtain a 
computationally more efficient method, we next discuss some practical imple- 
mentation issues. 

5.1 Sparsity patterns 

We start by discussing the sparsity pattern of the CDG method and compare 
it with that of the LDG and BR2 methods. We assume throughout that nodal 
bases [13] are used to span the approximating and weighting Galerkin spaces. 
For illustration purposes, we consider the triangular mesh shown in Figure 1, 
consisting of four elements and a finite element space of piecewise polynomials 
of degree p = 3 on each element. The total number of degrees of freedom is 
60, corresponding to 15 degrees of freedom per element. The sparsity patterns 
corresponding to the CDG, LDG, and BR2 methods are also shown in Figure 1. 
We note that the sparsity pattern of the IP method is identical to that of the 
BR2 method, and therefore the same remarks apply. 

As is well known, the LDG scheme introduces connections between degrees 
of freedom in nonneighboring elements. In this example, some degrees of free- 
dom in element 3 are connected to degrees of freedom in element 4. These 
connections are caused by the stabilization term (37), which involves the prod- 
uct of global Ufting functions. We note that these nonlocal connectivities also 
occur for quadrilateral discretizations and cannot be avoided by a more careful 
renumbering of the elements and/or internal interfaces [15]. 

Of the three schemes, the CDG method produces the smallest number of 
nonzero entries in the matrix. In fact, any nonzero entry in the CDG matrix is 
also a nonzero entry in the matrices generated by the other two schemes. The 
BR2 scheme is compact but connects the face nodes of each element with all 
the nodes of the neighboring element sharing that face. On the other hand, the 
CDG scheme connects only the nodes of those faces for which the switch is one, 
to the interior nodes of the neighboring element sharing that face. 
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Figure 1: The sparsity structure for four triangles with p = 3 (left plot). The 
CDG and the BR2 scheme are both compact in the sense that they connect 
only neighboring triangles; however, BR2 introduces more nonzeros. The LDG 
scheme is noncompact and gives connections between some nonneighboring tri- 
angles (3 and 4). 

5.2 Storage requirements 

In order to quantify the matrix storage requirements for the three schemes, we 
consider a simplex dement in d diinensions having d + 1 distinct neighboring 
elements. For polynomial basis function of degree p, the number of degrees 
of freedom per element is given hy S = (^^''), and the number of degrees of 
freedom along each element face is given by = (^^-1^)" Using this notation, 
we can obtain expressions for the number of nonzero matrix entries per interior 
element. 

For the CDG scheme we have one diagonal block with 5^ entries and d+1 
off-diagonal blocks with SeS entries. Since the scheme connects some element 
face nodes to all the nodes of the neighboring element sharing that face, we have 

McDG^S^ + {d+l)SeS. 

For the LDG scheme, the pattern is the same as for the CDG algorithm 
plus the additional nonlocal connectivities. Each such connectivity involves 
entries since the scheme connects face nodes to nonneighboring face nodes. The 
number of nonlocal connections a depends on the mesh and the switch, but 
on average, we have that in one dimension the switch can be chosen such that 
a = 0, and our experiments indicate that a « 1 for d = 2 and a i=a 2 for d = 3. 
The total number of nonzeros is then 

Mldg = S^ + {d+ l)SeS + aSl 
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Dim 






n — 9 
p — z, 


r) — 3 


p — 'i 


p — 5 


1 


CDG 


8 


15 


24 


35 


48 




LDG 


8 


15 


24 


35 


48 




BR2 


10 


19 


30 


43 


58 


2 


CDG 


27 


90 


220 


450 


819 




LDG 


31 


99 


236 


475 


855 




BR2 


33 


117 


292 


600 


1089 


3 


CDG 


64 


340 


1200 


3325 


7840 




LDG 


82 


412 


1400 


3775 


8722 




BR2 


76 


436 


1600 


4525 


10780 



Table 2: Memory requirements per interior simplex element for the CDG, LDG, 
and BR2 schemes. The case p = 3 in two dimensions is illustrated in Figure 1. 
The LDG scheme is assumed to have a = 0, 1, 2 noncompact neighbors in one, 
two, and three dimensions, respectively. 



Finally, for the BR2 (and also the IP) scheme, the pattern is the same as with 
the CDG scheme, but with the additional connections caused by the fact that 

all the face nodes connect to all the interior nodes in the neighboring elements. 
This results in S^S + Se{S — Se) entries per block, giving a total number of 
nonzeros of 

MBR2 = -S' + (d+l)(25-^e)5e. 

The memory requirements for d = 1,2,3 and p ~ 1, . . . , 5 are shown in Ta- 
ble 2. We note that the CDG method has the lowest memory requirements. 
For instance, in three dimensions with polynomials of degree p — A, the addi- 
tional storage requirements of the LDG and BR2 methods are 14% and 36%, 
respectively. 

Finally, we note that the CDG sparsity pattern is such that in addition 
to having fewer nonzero entries, the entire matrix can be stored using simple 
blockwise dense arrays. In particular, for a problem involving T elements, we 
can use an. S x S x T dense array for the diagonal blocks, and an S x Se x {d + 
1) X T dense array for the off-diagonal blocks. This representation is not only 
simple and compact, it also makes it straightforward to apply high-performance 
libraries such as the BLAS routines [4] for basic matrix operations. 

A similar storage format is harder to define for the LDG scheme, because 
of the noncompactness and the somewhat complex pattern in which these addi- 
tional blocks appear. For the BR2 scheme, while it is compact, and in principle 
one could use a storage scheme similar to that of the CDG method, the sparsity 
pattern of the off-diagonal blocks is nonrectangular, and therefore any dense 
storage strategy would require, at least, an additional array. 
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6 Numerical results 



In this section, we present some numerical experiments to assess the accuracy 
and behavior of the CDG algorithm. We consider a two-dimensional model 
problem. The domain is the unit square [0,1] x [0,1]. Dirichlet conditions 
are imposed at all the boundaries {dflo = dfl), and we choose the analytical 
solution 

u{x, y) = exp [asin(aa; + by) + (3cos{cx + dy)] (43) 

with numerical parameters a = 0.1, /3 = 0.3, a = 5.1, b = —6.2, c = 4.3, d = 3.4. 

We then solve the model Poisson problem (1) with the parameter k = 1 and 
with the Dirichlet boundary conditions gD{x,y) = u{x,y)\QQ,^. The source 
term, f{x,y), is obtained by analytical differentiation of (43). 

We consider triangular meshes obtained by splitting a regular nxn Cartesian 
grid into a total of 2n^ triangles, giving uniform element sizes oi h = 1/n. On 
these meshes, we consider solutions of polynomial degree p represented using a 
nodal basis within each triangle, with the nodes uniformly distributed. We use 
five different meshes, n = 2,4,8, 16,32, and five polynomial degrees, p = 1 to 
p = 5. 

6.1 Effect of the stabilization parameter Cn 

In order to assess the effect of the stabilization parameter, we discretize the 
Poisson equation (1) in two dimensions and solve for the numerical solution Uh 
using different values of the stabilization parameter Cn. The resulting equation 
system is solved using a preconditioned iterative solver [14]. We then compute 
the L2 error jju — ■u;i||o,n. The computed L2 error, \\u — Uh\\o,Q., is shown in 
Table 3 for the different values of p and n, and for Cn = 0, 1, and 10, using 
the consistent switch. The same results are reported for the natural switch in 
Table 4. 

Wc note that the accuracy is only weakly dependent on the value of Cn. 
The only noticeable differences are for the underresolved cases {p = 1,2 and 
n = 2) when using a large amount of stabilization, Cn = 10. We obtain the 
optimal convergence rate of p-l- 1 for all cases. Using the natural switch, instead 
of the consistent one, makes the errors somewhat larger, but on average only by 
11% and, in the worst case, only by 42%. 

Table 5 shows the errors and the convergence rates for the gradient of the 
solution using the CDG method with Cn = 0. In particular, we calculate 
the seminorm {^xeT^ 1^ ~ )^^'^- observe optimal convergence at the 

expected rate of p. 

6.2 Comparison with the LDG and BR2 schemes 

Here, we discretize the equations using the CDG, LDG, and BR2 schemes. For 
the CDG and the LDG methods, we use the consistent switch and set Cn = 0, 
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V 




n 


= 2 




n 


= 4 




n 


= 8 




n = 


= 16 




n = 


= 32 




Rate 


1 





4.55 


• 10- 


2 


1.52 


10- 


-2 


4.63 


10 


-3 


1.26- 


10- 


3 


3.27 


10- 


4 


1.9 




1 


4.55 


• 10- 


2 


1.49 


10- 


-2 


4.56 


10 


-3 


1.25- 


10- 


3 


3.26 


10- 


4 


1.9 




10 


2.20 


• 10- 





2.07 


10- 


-2 


4.24 


10 


-3 


1.16- 


10- 


3 


3.13 


10- 


4 


1.9 


2 





9.00 


• 10- 


3 


1.80 


10- 


-3 


2.56 


10 


-4 


3.36- 


10- 


5 


4.29 


10- 


6 


3.0 




1 


9.10 


• 10- 


3 


1.80 


10- 


-3 


2.56 


10 


-4 


3.36 • 


10- 


5 


4.29 


10- 


6 


3.0 




10 


2.89 


• 10- 


2 


2.01 


10- 


-3 


2.62 


10 


-4 


3.38- 


10- 


5 


4.30 


10- 


6 


3.0 


3 





2.61 


• 10- 


3 


2.44 


10- 


-4 


1.72 


10 


-5 


1.11 • 


10- 


6 


7.04 


10- 


8 


4.0 




1 


2.63 


• 10- 


3 


2.44 


10- 


-4 


1.72 


10 


-5 


1.11 • 


10- 


6 


7.04 


10- 


8 


4.0 




10 


4.16 


• 10- 


3 


2.59 


10- 


-4 


1.73 


10 


-5 


1.11 • 


10- 


6 


7.03 


10- 


8 


4.0 


4 





1.09 


• 10- 


3 


4.52 


10- 


-5 


1.57 


10 


-6 


5.14- 


10- 


8 


1.64 


10- 


9 


5.0 




1 


1.09 


• 10- 


3 


4.54 


10- 


-5 


1.57 


10 


-6 


5.15- 


10- 


8 


1.64 


10- 


9 


5.0 




10 


1.19 


• 10- 


3 


4.77 


10- 


-5 


1.60 


10 


-6 


5.16- 


10- 


8 


1.64 


10- 


9 


5.0 


5 





3.73 


• 10- 


4 


9.31 


10- 


-6 


1.76 


10 


-7 


2.83- 


10- 


9 


4.47 


10- 


11 


6.0 




1 


3.75 


• 10- 


4 


9.32 


10- 


-6 


1.76 


10 


-7 


2.83 • 


10- 


9 


4.47 


10- 


11 


6.0 




10 


4.07 


• 10- 


4 


9.52 


10- 


-6 


1.77 


10 


-7 


2.84- 


10- 


9 


4.47 


10- 


11 


6.0 



Table 3: L2 errors in the solution for the model Poisson problem, for various 
polynomial degrees p, mesh sizes n, and Cn values. The consistent switch is 
used. The convergence rate is calculated based on the two finest meshes. 



19 



V 


Cii 


n 


= 2 




n 


= 4 




n 


= 8 




n = 


= 16 




n 


= 32 




Rate 


1 





3.72 


10" 


2 


1.61 


• 10- 


2 


4.71 


• 10 


-3 


1.30 • 


10- 


3 


3.39 


• 10- 


4 


1.9 




1 


3.83 


10" 


2 


1.50 


• 10- 


2 


4.70 


• 10 


-3 


1.32- 


10- 


3 


3.38 


• 10- 


4 


2.0 




10 


2.33 


10- 


1 


3.40 


• 10- 


2 


4.64 


• 10 


-3 


1.25- 


10- 


3 


3.31 


• 10- 


4 


1.9 


2 





1.28 


10" 


2 


1.96 


• 10- 


3 


3.03 


• 10 


-4 


3.98 • 


10- 


5 


5.04 


• 10- 


6 


3.0 




1 


1.18 


10- 


2 


2.07 


• 10- 


3 


2.88 


• 10 


-4 


4.01 • 


10- 


5 


5.02 


• 10- 


6 


3.0 




10 


3.24 


10- 


2 


3.00 


• 10- 


3 


3.37 


• 10 


-4 


4.05- 


10- 


5 


5.16 


• 10- 


6 


3.0 


3 





3.03 


10" 


3 


2.68 


• 10- 


4 


2.01 


• 10 


-5 


1.33 • 


10- 


6 


8.63 


• 10- 


8 


4.0 




1 


3.25 


10- 


3 


2.74 


• 10- 


4 


2.05 


• 10 


-5 


1.33- 


10- 


6 


8.61 


• 10- 


8 


3.9 




10 


1.84 


10- 


2 


3.56 


• 10- 


4 


2.28 


• 10 


-5 


1.38- 


10- 


6 


8.79 


• 10- 


8 


4.0 


4 





9.67 


10- 


4 


5.15 


• 10- 


5 


1.82 


• 10 


-6 


5.86- 


10- 


8 


1.87 


• 10- 


9 


5.0 




1 


1.33 


10- 


3 


5.24 


• 10- 


5 


1.81 


• 10 


-6 


5.90- 


10- 


8 


1.88 


• 10- 


9 


5.0 




10 


1.98 


10- 


3 


6.30 


• 10- 


5 


1.95 


• 10 


-6 


6.12- 


10- 


8 


1.90 


• 10- 


9 


5.0 


5 





3.98 


10- 


4 


1.01 


• 10- 


5 


1.85 


• 10 


-7 


3.07- 


10- 


9 


4.83 


10- 


11 


6.0 




1 


3.84 


10- 


4 


1.02 


• 10- 


5 


1.88 


• 10 


-7 


3.06 • 


10- 


9 


4.85 


10- 


11 


6.0 




10 


4.97 


10- 


4 


1.15 


• 10- 


5 


1.98 


• 10 


-7 


3.11 • 


10- 


9 


4.88 


10- 


11 


6.0 



Table 4: L2 errors in the solution for the model Poisson problem, for various 
polynomial degrees p, mesh sizes n, and Cn values. The natural switch is used. 



p 


n 


= 2 




n 


= 4 




n 


= 8 




n = 


16 




n 


= 32 


Rate 


1 


1.80 


• 10- 


-0 


6.09 


• 10- 


1 


3.05 


• 10- 


-1 


1.54- 


10- 


-1 


7.75 


• 10-2 


1.0 


2 


7.40 


• 10- 


-1 


1.57 


• 10- 


1 


3.73 


• 10- 


-2 


9.20 • 


10- 


-3 


2.28 


• 10-^ 


2.0 


3 


2.57 


• 10- 


-1 


3.01 


• 10- 


2 


3.63 


• 10- 


-3 


4.37- 


10- 


-4 


5.36 


• 10-'^ 


3.0 


4 


9.53 


• 10- 


-2 


5.96 


• 10- 


3 


3.61 


• 10- 


-4 


2.18 • 


10- 


-5 


1.32 


• 10-6 


4.0 


5 


5.42 


• 10- 


-2 


1.33 


• 10- 


3 


3.67 


• 10- 


-5 


1.04 • 


10- 


-6 


3.11 


• 10-8 


5.0 



Table 5: The errors in the gradient for the CDG scheme with consistent switch 
and Cii = 0. 



20 




except at the Dirichlct boundaries, where C'u = 1. The hfting parameter in the 
BR2 scheme is 77 = 3, which is the value required for stabihty [5]. 

The accuracy results for the CDG, LDG, and BR2 schemes are shown in 
Figure 2, with details in Table 6. We note that the CDG Hc;li{^ni{^ is the most 
accurate scheme in most of the test cases. For low polynomial degrees and on 
the coarse meshes, the difference is often more than a factor of 2, while for 
well-resolved solutions, CDG and LDG are similar, and BR2 is about 10% less 
accurate. We can also see that all schemes give optimal convergence rates close 
to p+ 1 for it/i||o,n- 

6.3 Spectral radius 

In our next study, we compute the spectral radius |Aniax| of the discretized 
matrix and compare the three methods. The spectral radius of the discretized 
matrix determines the magnitude of the timestep if an explicit time marching 
solution is sought. In Table 7, we show these values for each of the simulations 
in the previous section, scaled by the factor (h/p)^. Here we have used the 
consistent switch with the constant Cn = for the CDG and LDG methods 
and a value of ry = 3 in the BR2 discretization. We observe that the CDG and the 
LDG methods have almost identical spectral radii, while the BR2 method gives 
about 50% larger values. It is possible that a lower value of the rj parameter in 
the BR2 method may reduce the spectral radius. However, in this case stability 
may be compromised. 
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p 


Scheme 


n 


= 2 




n 


= 4 




n 


= 8 




n = 


= 16 




n 


= 32 




Rate 


1 


CDG 


4.54 


• 10- 


2 


1.52 


• 10- 


2 


4.62 


• 10- 


3 


1.25 • 


10- 


3 


3.27 


• 10- 


4 


1.9 




LDG 


1.34 


• 10- 


1 


1.73 


• 10- 


2 


4.68 


• 10- 


3 


1.25 • 


10- 


3 


3.26 


• 10- 


4 


1.9 




BR2 


8.60 


• 10- 


2 


3.08 


• 10- 


2 


9.23 


• 10- 


3 


2.47- 


10- 


3 


6.36 


• 10- 


4 


2.0 


2 


CDG 


8.99 


• 10- 


3 


1.79 


• 10- 


3 


2.55 


• 10- 


4 


3.35- 


10- 


5 


4.28 


• 10- 


6 


3.0 




LDG 


3.81 


• 10- 


2 


2.92 


• 10- 


3 


3.03 


• 10- 


4 


3.59- 


10- 


5 


4.42 


• 10- 


6 


3.0 




BR2 


1.66 


• 10- 


2 


2.75 


• 10- 


3 


3.16 


• 10- 


4 


3.75- 


10- 


5 


4.60 


• 10- 


6 


3.0 


3 


CDG 


2.61 


• 10- 


3 


2.44 


• 10- 


4 


1.71 


• 10- 


5 


1.10- 


10- 


6 


7.03 


• 10- 


8 


4.0 




LDG 


5.88 


• 10- 


3 


3.81 


• 10- 


4 


2.04 


• 10- 


5 


1.18- 


10- 


6 


7.23 


• 10- 


8 


4.0 




BR2 


5.64 


• 10- 


3 


3.77 


• 10- 


4 


2.47 


• 10- 


5 


1.52- 


10- 


6 


9.46 


• 10- 


8 


4.0 


4 


CDG 


1.09 


• 10- 


3 


4.52 


• 10- 


5 


1.56 


• 10- 


6 


5.14- 


10- 


8 


1.63 


• 10- 


9 


5.0 




LDG 


2.04 


• 10- 


3 


5.00 


• 10- 


5 


1.65 


• 10- 


6 


5.28- 


10- 


8 


1.66 


• 10- 


9 


5.0 




BR2 


1.30 


• 10- 


3 


6.22 


• 10- 


5 


2.05 


• 10- 


6 


6.57- 


10- 


8 


2.07 


• 10- 


9 


5.0 


5 


CDG 


3.73 


• 10- 


4 


9.30 


• 10- 


6 


1.75 


• 10- 


7 


2.83- 


10- 


9 


4.46 


10- 


11 


6.0 




LDG 


1.06 


• 10- 


3 


1.32 


• 10- 


5 


1.93 


• 10- 


7 


2.91 • 


10- 


9 


4.50 


10- 


11 


6.0 




BR2 


4.42 


• 10- 


4 


1.08 


• 10- 


5 


2.05 


• 10- 


7 


3.31 • 


10- 


9 


5.23 


10- 


11 


6.0 



Table 6: L2 errors in the solution for the model Poisson problem, for different 
polynomial degree, p, and mesh size, n, using CDG, LDG, and BR2 Schemes. 
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4 
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;( = 32 


1 


CDG 


153 


4 


157 


5 


159 


4 


159.9 


160.1 




LDG 


149 


5 


156 


7 


159 


2 


159.9 


160.1 




BR2 


244 





244 


8 


245 


2 


245.4 


245.4 


2 


CDG 


137 


4 


139 


8 


140 


8 


141.1 


141.1 




LDG 


135 


1 


139 


5 


140 


7 


141.1 


141.1 




BR2 


216 


1 


215 


5 


215 


3 


215.1 


215.1 


3 


CDG 


159 


9 


161 


3 


161 


8 


162.0 


162.0 




LDG 


159 


5 


161 


1 


161 


8 


162.0 


162.0 




BR2 


244 


4 


244 





243 


8 


243.8 


243.8 


4 


CDG 


198 


4 


200 


3 


201 





201.2 


201.3 




LDG 


197 


7 


200 


2 


201 





201.2 


201.3 




BR2 


302 


1 


300 


9 


300 


6 


300.6 


300.6 


5 


CDG 


244 


8 


246 





246 


4 


246.5 


246.5 




LDG 


245 


1 


246 





246 


4 


246.5 


246.5 




BR2 


368 


5 


368 


4 


368 


4 


368.4 


368.4 



Table 7: The spectral radii of the matrices for the model Poisson problem, scaled 
by {h/pf. 
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7 Conclusions 



We have presented a new scheme for discretizing elliptic operators in the context 
of discontinuous Galerkin approximations. The main advantage of the proposed 
scheme is its reduced sparsity pattern when compared to alternative schemes 
such as the LDG, BR2, or IP methods. This is important when an implicit solu- 
tion technique is required. Compared to the LDG scheme the proposed scheme 
is compact, meaning that only degrees of freedom in neighboring elements are 
connected. Compared to the BR2 and IP schemes, which are also compact, the 
CDG scheme produces a smaller number of nonzero entries in the off-diagonal 
blocks and, at the same time, the nonzero elements in the CDG scheme are 
amenable to a dense block matrix storage. Like the alternative approaches, 
the proposed scheme converges optimally, and numerical tests indicate that the 
accuracy obtained compares well with that of the LDG or BR2 schemes. An 
additional potential advantage of the CDG scheme over the LDG scheme when 
both schemes are used with minimal dissipation (i.e., Cn = in the interior 
faces) is its insensitivity to the face ordering. 
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